clear all
set more off
cap log close

********************************************************************************
***** Project: The Short and Long Term Effects of In-Person Performance Feedback
********************************************************************************
***** A. R. Soetevent & G. J. Romensen
********************************************************************************
***** Basic regressions - Coaching: Heterogeneous Treatment Response
********************************************************************************

*global filepath "C:\JPEMicReplication"
*global paperpath "$filepath\TablesGraphs"

local abcd "acceleratie rem bochten fueleconomyLpKM"
log using "$filepath/Logs/X08HetEffectCoachingProductivity.log", replace

/*** Notes ***/
*=> No fuel economy observations postfeedback in urban area.
*=> Fuel Economy: kilometers per liter of fuel
*=> ABC dimensions: number of events per 10 kilometers
/*************/

use "$filepath\DEPO\DataMainAnalysisDEPO.dta"

*** Merge with data on EB-adjusted driver fixed effects

sort chauf_nr_rug
merge m:1 chauf_nr_rug using "$filepath\DEPO\EstimatedFixedEffectsDriverLevelDEPO.dta", keepusing(feEBacceleratie feEBrem feEBbochten feEBfueleconomyLpKM feEBrankacceleratie feEBrankrem feEBrankbochten feEBrankfueleconomyLpKM IsCoach)
drop _merge

rename bustype bustypes
** Drop all Irisbus observations 
drop if bustypes==3
gen byte ZH = 0 
replace ZH = 1 if regio == "ZH"

** Analysis: based on Treatment region
drop if regio == "ZH"

replace intouro = 0 if intouro ==.

* Drop months with imperfect tracking by coaches 
drop if datum>date("30-4-2016", "DMY")

* Drop eco-coaches [chauf_nr_rug: randomly generated depository numbers!]
local eco_nr "939 1404 519 1286 1610 531"

foreach x of local eco_nr {
drop if chauf_nr_rug==`x'
}


** Newly defined variables [geplande_ritafstand is in meters, for this reason *1000]
gen aantal_haltespKM = (1000*aantal_haltes)/geplande_ritafstand
label variable aantal_haltespKM "Nr. stops per km."


****************************************
*** A. Create global list of covariates
****************************************
* outcome variables
global abcd "acceleratie rem bochten fueleconomyLpKM"

global covMaand "maand_1 maand_2 maand_3 maand_4 maand_6 maand_7 maand_8 maand_9 maand_10 maand_11 maand_12" 
* (May = default)

global covBusType "vdl10 vdl14  iris10 iris10cng iris12 iris12cng intouro"
* (vdl12 = default - more than 50 per cent of observations)

global covWeatherExt "Temp6_3C Temp3_0C Temp0_3C Temp3_6C Temp9_12C Temp12_15C Temp15_18C Temp18_21C Temp21_24C Temp24_27C Rain5_10mm Rain10_15mm Rain15_20mm Rain20_30mm Rain30_40mm Wind0_2ms Wind4_6ms Wind6_8ms Wind8_10ms Wind10_20ms"  
*(Defaults: Wind2_4ms, Temp6_9C, Rain0_5mm)
global covWeather "Temp10_5C Temp5_10C Temp15_20C Temp20_25 Rain0_5mm Rain5_50mm Wind0_3ms Wind5_8ms Wind8_20ms"  
*(Defaults: Wind4_5ms, Temp10_15C, Rain0_0mm)


gen triplength=geplande_ritafstand/1000
global covEnvironm "ochtendspits avondspits uitleenrit aantal_haltespKM stadsrit triplength"
** Note: this is other covEnvironm than used in other do files: "geplande_ritafstand aantal_haltes" have been replaced by "aantal_haltespKM"

global covPassengers "lnovcheckins ovcheckinsmissing"

global covEndogenous "punctuality"


** Postcoaching dummies: months (21 day periods)
gen byte MonthsToCoachingplus0 = 0
foreach x of numlist 1(1)3 {
	gen byte MonthsToCoachingmin`x' = 0
	gen byte MonthsToCoachingplus`x' = 0
	gen byte LOWMonthsToCoachingmin`x' = 0
	gen byte LOWMonthsToCoachingplus`x' = 0
	gen byte MIDMonthsToCoachingmin`x' = 0
	gen byte MIDMonthsToCoachingplus`x' = 0
	gen byte HIGHMonthsToCoachingmin`x' = 0
	gen byte HIGHMonthsToCoachingplus`x' = 0
}

gen byte MonthsToCoachingminbin = 0 
gen byte MonthsToCoachingplusbin = 0 

replace  MonthsToCoachingminbin = 1 if  (daysbeforecoachingA70 == 1 | daysbeforecoachingA70min == 1)
replace  MonthsToCoachingmin3 = 1 if    (daysbeforecoachingA63 == 1 | daysbeforecoachingA56 == 1 | daysbeforecoachingA49 == 1)
replace  MonthsToCoachingmin2 = 1 if    (daysbeforecoachingA42 == 1 | daysbeforecoachingA35 == 1 | daysbeforecoachingA28 == 1)
replace  MonthsToCoachingmin1  = 1 if  (daysbeforecoachingA7 == 1 | daysbeforecoachingA14 == 1 | daysbeforecoachingA21 == 1)
replace  MonthsToCoachingplus0  = 1 if  (daysaftercoachingA0 == 1)
replace  MonthsToCoachingplus1 = 1 if  (daysaftercoachingA7 == 1 | daysaftercoachingA14 == 1 | daysaftercoachingA21 == 1)
replace  MonthsToCoachingplus2 = 1 if  (daysaftercoachingA28 == 1 | daysaftercoachingA35 == 1 | daysaftercoachingA42 == 1)
replace  MonthsToCoachingplus3 = 1 if  (daysaftercoachingA49 == 1 | daysaftercoachingA56 == 1 | daysaftercoachingA63 == 1)
replace  MonthsToCoachingplusbin  = 1 if  (daysaftercoachingA70 == 1 | daysaftercoachingA70plus == 1)

gen byte MonthsToCoaching = . 
replace MonthsToCoaching  = - 3  if  MonthsToCoachingmin3  == 1
replace MonthsToCoaching  = - 2  if  MonthsToCoachingmin2  == 1
replace MonthsToCoaching  = - 1  if  MonthsToCoachingmin1  == 1
replace MonthsToCoaching  =   0  if  MonthsToCoachingplus0 == 1
replace MonthsToCoaching  =   1  if  MonthsToCoachingplus1 == 1
replace MonthsToCoaching  =   2  if  MonthsToCoachingplus2 == 1
replace MonthsToCoaching  =   3  if  MonthsToCoachingplus3 == 1

** Does driver have trip observations in all of MonthsToCoaching  = - 1, MonthsToCoaching  = 0, MonthsToCoaching  = + 1, MonthsToCoaching  = + 2, MonthsToCoaching  = + 3? 


foreach var in $abcd {
		gen byte LOWMin_`var' = (feEBrank`var' >=75)
		gen byte MIDMin_`var' = ((feEBrank`var' < 75) & (feEBrank`var' > 25))
		gen byte HIGHMin_`var' = (feEBrank`var' <= 25)
	
	foreach x of numlist 3(1)1 {
		gen byte LOWMin`x'_`var' = MonthsToCoachingmin`x' * LOWMin_`var'
		gen byte MIDMin`x'_`var' = MonthsToCoachingmin`x' * MIDMin_`var'
		gen byte HIGHMin`x'_`var' = MonthsToCoachingmin`x' * HIGHMin_`var'
	}
	foreach x of numlist 0(1)3 {
		gen byte LOWPlus`x'_`var' = MonthsToCoachingplus`x' * (feEBrank`var' >=75)
		gen byte MIDPlus`x'_`var' = MonthsToCoachingplus`x' * ((feEBrank`var' < 75) & (feEBrank`var' > 25))
		gen byte HIGHPlus`x'_`var' = MonthsToCoachingplus`x' * (feEBrank`var' <= 25)
	}

		gen byte LOWMinbin_`var' = MonthsToCoachingminbin * (feEBrank`var' >=75)
		gen byte MIDMinbin_`var' = MonthsToCoachingminbin * ((feEBrank`var' < 75) & (feEBrank`var' > 25))
		gen byte HIGHMinbin_`var' = MonthsToCoachingminbin * (feEBrank`var' <= 25)

		gen byte LOWPlusbin_`var' = MonthsToCoachingplusbin * (feEBrank`var' >=75)
		gen byte MIDPlusbin_`var' = MonthsToCoachingplusbin * ((feEBrank`var' < 75) & (feEBrank`var' > 25))
		gen byte HIGHPlusbin_`var' = MonthsToCoachingplusbin * (feEBrank`var' <= 25)
}


** coaching dummies
* Post
global  covPostCoachingMonth
	foreach x of numlist 0(1)3 {
		global covPostCoachingMonth "$covPostCoachingMonth MonthToCoachingplus`x'" 
	}
	global covPostCoachingMonth "$covPostCoachingMonth MonthToCoachingplusbin" 
** Pre[month before coaching is baseline]
global covPreCoachingMonth
global covPreCoachingMonth "$covPreCoachingMonth MonthsToCoachingminbin"
	foreach x of numlist 3(1)2 {
		global covPreCoaching "$covPreCoaching MonthsToCoachingmin`x'" 
	}
	
* Het. coaching dummies	
foreach var in $abcd {
	global  covCoachingMonth_`var' "LOWMinbin_`var'  LOWMin3_`var' LOWMin2_`var'  LOWPlus0_`var'  LOWPlus1_`var'  LOWPlus2_`var' LOWPlus3_`var'  LOWPlusbin_`var' MIDMinbin_`var'  MIDMin3_`var' MIDMin2_`var'  MIDPlus0_`var'  MIDPlus1_`var'  MIDPlus2_`var' MIDPlus3_`var'  MIDPlusbin_`var' HIGHMinbin_`var' HIGHMin3_`var'  HIGHMin2_`var'  HIGHPlus0_`var'  HIGHPlus1_`var'  HIGHPlus2_`var'  HIGHPlus3_`var'  HIGHPlusbin_`var'"
	global  covMonth  "MonthsToCoachingminbin  MonthsToCoachingmin3 MonthsToCoachingmin2 MonthsToCoachingplus0 MonthsToCoachingplus1 MonthsToCoachingplus2 MonthsToCoachingplus3 MonthsToCoachingplusbin"
}



** Postcoaching dummies: non treatment-specific
global  covPostCoaching
	foreach x of numlist 0(7)70 {
		global covPostCoaching "$covPostCoaching daysaftercoachingA`x'" 
	}
	global covPostCoaching "$covPostCoaching daysaftercoachingA70plus"

	** Postcoaching dummies: treatment-specific
global  covPostCoachingT
foreach i of numlist 1(1)4 {
	foreach x of numlist 0(7)70 {
		global covPostCoachingT "$covPostCoachingT daysaftercoaching`x'T`i'" 
	}
	global covPostCoachingT "$covPostCoachingT daysaftercoaching70T`i'plus"
}

** Precoaching dummies: non treatment-specific
global covPreCoaching
global covPreCoaching "$covPreCoaching daysbeforecoachingA70min"
	foreach x of numlist 70(7)14 {
		global covPreCoaching "$covPreCoaching daysbeforecoachingA`x'" 
	}

********************************************************************************
*** B. Time trend in ABC + Fuel economy correcting for time-variant variables
********************************************************************************
set matsize 10000

 
/*** Determine the set of observations used for the analysis ***/
gen byte regobsfuel=1
replace regobsfuel=0 if geplande_ritafstand==. | lnovcheckins==. | punctuality==. | aantal_haltespKM==. | dep_fueleconomyLpKM==.


gen byte regobsabc=1
replace regobsabc=0  if geplande_ritafstand==. | lnovcheckins==. | punctuality==. | aantal_haltespKM==. | dep_acceleratie==. | dep_bochten==. | dep_rem==.
xtset chauf_nr_rug 

** No. of observation per driver per outcome dimension pre-coaching
sort chauf_nr_rug 
by chauf_nr_rug: egen postcoachingobs = total(postcoaching == 1)
gen regobsabcprecoaching = regobsabc * (1-postcoaching)
by chauf_nr_rug: egen NobsPerDriverdep_abc = total(regobsabcprecoaching == 1)
gen regobsfuelprecoaching = regobsfuel * (1-postcoaching)
by chauf_nr_rug: egen NobsPerDriverdep_fuel = total(regobsfuelprecoaching == 1)
drop regobsabcprecoaching regobsfuelprecoaching

** Summary statistics het. treatment effects
foreach var in $abcd {
	tabstat dep_`var'  if LOWMin_`var' == 1, s(mean semean n) by(MonthsToCoaching)
    tabstat dep_`var'  if MIDMin_`var' == 1, s(mean semean n) by(MonthsToCoaching)
	tabstat dep_`var'  if HIGHMin_`var' == 1, s(mean semean n)by(MonthsToCoaching)
}

foreach var in $abcd {
preserve
	collapse (mean) dep_`var', by(chauf_nr_rug MonthsToCoaching LOWMin_`var' MIDMin_`var' HIGHMin_`var') cw
	drop if MonthsToCoaching ==. 
	bysort chauf_nr_rug: egen MtCmin2Present = max(MonthsToCoaching == -2)
	bysort chauf_nr_rug: egen MtCmin1Present = max(MonthsToCoaching == -1)
	bysort chauf_nr_rug: egen MtC0Present = max(abs(MonthsToCoaching) <= 0.01)
	bysort chauf_nr_rug: egen MtC1Present = max(MonthsToCoaching == 1)
	bysort chauf_nr_rug: egen MtC2Present = max(MonthsToCoaching == 2)
	
	* only include drivers with observations (=trips) in month prior to coaching, month of coaching + the two months after.
	gen completeset = (MtCmin1Present + MtC0Present + MtC1Present + MtC2Present == 4)
	drop if completeset == 0
	
	collapse (mean) mean_`var' = dep_`var' (semean) semean_`var' = dep_`var' (count) n_`var' = dep_`var' , by(MonthsToCoaching LOWMin_`var' MIDMin_`var' HIGHMin_`var') 
	sort LOWMin_`var' MIDMin_`var' HIGHMin_`var' MonthsToCoaching 
	save "$paperpath\tempMtC_`var'.dta", replace
	
	** Graph
	gen lowCI = mean_`var' - 1.96*semean_`var'
	gen upCI = mean_`var' + 1.96*semean_`var'

	drop if MonthsToCoaching < -1 |  MonthsToCoaching > 3
twoway 	(scatter mean_`var' MonthsToCoaching if LOWMin_`var' == 1  , ms(Oh) mc(red)) (rcap lowCI upCI MonthsToCoaching if LOWMin_`var' == 1 , lc(red))   ///
		(scatter mean_`var' MonthsToCoaching if MIDMin_`var' == 1  , ms(Oh) mc(blue)) (rcap lowCI upCI MonthsToCoaching if MIDMin_`var' == 1 , lc(blue)) ///
		(scatter mean_`var' MonthsToCoaching if HIGHMin_`var' == 1  , ms(Oh) mc(green)) (rcap lowCI upCI MonthsToCoaching if HIGHMin_`var' == 1 , lc(green)) ///
       , name(p`var', replace) 	   ///
	   legend(lab(1 "Lower 25%") lab(2 "") lab(3 "Middle 50%") lab(4 "") lab(5 "Upper 25%") lab(6 "") )  ///
	   xscale(r(-1.1 3.1))  ///
	   xtitle("Months to Coaching") ytitle("`var'")   ///
	   title(" ", size(medium))  ///
	   xlab(-1 "-1" 0 "0" 1 "1" 2 "2" 3 "3" )
	
restore
}


** REGRESSIONS **
gen byte selectie=0
local i = 0
foreach var in $abcd {
	local i = `i' + 1
		if `i'<4 {
		replace selectie=regobsabc if  NobsPerDriverdep_abc >= 25
		}
		if `i'==4 {
		replace selectie=regobsfuel if NobsPerDriverdep_fuel >= 25
			}
		
		di "Het. Fixed effects regression: `var'"
		
		** LOWEST 25%
		* (1) FE - without controls + day and bustype fixed effects (including interactions)
		xtreg dep_`var' postcoaching  i.datum if selectie==1 & LOWMin_`var' == 1, fe cluster(chauf_nr_rug)
		estadd loc controls "No", replace
		estadd loc driverfe "No", replace
		estadd loc dayfe 	"Yes", replace	
		estadd loc busdayfe "No", replace 
		estimates store hteLOW`var'1,  title(hteLOW`var'1)
		* (2) FE + day and bustype fixed effects (including interactions) - Non treatment specific impact coaching
		xtreg dep_`var' $covMonth  $covBusType $covPassengers $covEnvironm  i.bustypes i.datum i.bustypes#i.datum if selectie==1 & LOWMin_`var' == 1, fe cluster(chauf_nr_rug)
		*xtreg dep_`var' $covPreCoaching $covPostCoaching  $covBusType $covPassengers $covEnvironm  i.bustypes i.datum i.bustypes#i.datum if selectie==1 & LOWMin_`var' == 1, fe cluster(chauf_nr_rug)	
		estadd loc controls "Yes", replace
		estadd loc driverfe "Yes", replace
		estadd loc dayfe 	"Yes", replace	
		estadd loc busdayfe "Yes", replace 
		estimates store hteLOW`var'2,  title(hteLOW`var'2)
		
		** MIDDLE 50%
		* (1) FE - without controls + day and bustype fixed effects (including interactions)
		xtreg dep_`var' postcoaching  i.datum if selectie==1 & MIDMin_`var' == 1, fe cluster(chauf_nr_rug)
		estadd loc controls "No", replace
		estadd loc driverfe "No", replace
		estadd loc dayfe 	"Yes", replace	
		estadd loc busdayfe "No", replace 
		estimates store hteMID`var'1,  title(hteMID`var'1)
		* (2) FE + day and bustype fixed effects (including interactions) - Non treatment specific impact coaching
		xtreg dep_`var' $covMonth  $covBusType $covPassengers $covEnvironm  i.bustypes i.datum i.bustypes#i.datum if selectie==1 & MIDMin_`var' == 1, fe cluster(chauf_nr_rug)
		*xtreg dep_`var' $covPreCoaching $covPostCoaching  $covBusType $covPassengers $covEnvironm  i.bustypes i.datum i.bustypes#i.datum if selectie==1 & MIDMin_`var' == 1, fe cluster(chauf_nr_rug)	
		estadd loc controls "Yes", replace
		estadd loc driverfe "Yes", replace
		estadd loc dayfe 	"Yes", replace	
		estadd loc busdayfe "Yes", replace 
		estimates store hteMID`var'2,  title(hteMID`var'2)
		
		** HIGHEST 25%
		* (1) FE - without controls + day and bustype fixed effects (including interactions)
		xtreg dep_`var' postcoaching i.datum if selectie==1 & HIGHMin_`var' == 1, fe cluster(chauf_nr_rug)
		estadd loc controls "No", replace
		estadd loc driverfe "No", replace
		estadd loc dayfe 	"Yes", replace	
		estadd loc busdayfe "No", replace 
		estimates store hteHIGH`var'1,  title(hteHIGH`var'1)
		* (2) FE + day and bustype fixed effects (including interactions) - Non treatment specific impact coaching
		xtreg dep_`var'  $covMonth       $covBusType $covPassengers $covEnvironm  i.bustypes i.datum i.bustypes#i.datum if selectie==1 & HIGHMin_`var' == 1, fe cluster(chauf_nr_rug)
		*xtreg dep_`var' $covPreCoaching $covPostCoaching  $covBusType $covPassengers $covEnvironm  i.bustypes i.datum i.bustypes#i.datum if selectie==1 & HIGHMin_`var' == 1, fe cluster(chauf_nr_rug)	
		estadd loc controls "Yes", replace
		estadd loc driverfe "Yes", replace
		estadd loc dayfe 	"Yes", replace	
		estadd loc busdayfe "Yes", replace 
		estimates store hteHIGH`var'2,  title(hteHIGH`var'2)
}


***********************************************************************************
***  [Figure K.11: Heterogeneous Treatment Effects of Coaching - Treatment Region]
***********************************************************************************

local i = 0
foreach var in $abcd {
	local i = `i' + 1
	if `i'==1 { 
		local ynaam  "Change in number of events per 10km" 
		local yschaal "r(-2 1)" 
		local ylabeltje "-2 -1.5 -1 -0.5 0 0.5 1,labsize(*1.4)"
	}
	else if `i'==2 { 
			local ynaam  "Change in number of events per 10km" 
		local yschaal "r(-0.4 0.2)" 
		local ylabeltje "-0.4 -0.3 -0.2 -0.1 0 0.1 0.2, labsize(*1.4)"
	}
	else if `i'==3 { 
		local ynaam  "Change in number of events per 10km" 
		local yschaal "r(-1.0 0.2)" 
		local ylabeltje "-1.0 -0.8 -0.6 -0.4 -0.2 0 0.2, labsize(*1.4)"
	}
	else if `i'==4 { 
		local ynaam  "Change in fuel economy (liters/100km)" 
		local yschaal "r(-1.6 0.4)" 
		local ylabeltje "-1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0 0.2, labsize(*1.4)"
	}
	

coefplot (hteLOW`var'2), keep(MonthsToCoachingmin2 MonthsToCoachingplus0 MonthsToCoachingplus1 MonthsToCoachingplus2 MonthsToCoachingplus3) base vertical yline(0) xlabel("") xtick(#5)  xlabel(1 "-3 to -1" 2 "0" 3 "1 to 3" 4 "4 to 6" 5 "7 to 9" , labsize(*1.4)) ///
yscale(`yschaal') ylabel(`ylabeltje') ///
xtitle("weeks relative to coaching",  size(*1.4)) ytitle("`ynaam'" ,  size(*1.4)) graphregion(fcolor(white) lcolor(white)) bgcolor(white) ///
name(hteLOW_`var', replace)

coefplot (hteMID`var'2), keep(MonthsToCoachingmin2 MonthsToCoachingplus0 MonthsToCoachingplus1 MonthsToCoachingplus2 MonthsToCoachingplus3) vertical yline(0) xlabel("") xtick(#5)  xlabel(1 "-3 to -1" 2 "0" 3 "1 to 3" 4 "4 to 6" 5 "7 to 9" , labsize(*1.4)) ///
yscale(`yschaal') ylabel(`ylabeltje') ///
xtitle("weeks relative to coaching",  size(*1.4)) ytitle("`ynaam'",  size(*1.4)) graphregion(fcolor(white) lcolor(white)) bgcolor(white) ///
name(hteMID_`var', replace)

coefplot (hteHIGH`var'2), keep(MonthsToCoachingmin2 MonthsToCoachingplus0 MonthsToCoachingplus1 MonthsToCoachingplus2 MonthsToCoachingplus3) vertical yline(0) xlabel("") xtick(#5)  xlabel(1 "-3 to -1" 2 "0" 3 "1 to 3" 4 "4 to 6" 5 "7 to 9",  labsize(*1.4)) ///
yscale(`yschaal') ylabel(`ylabeltje') ///
xtitle("weeks relative to coaching",  size(*1.4)) ytitle("`ynaam'",  size(*1.4)) graphregion(fcolor(white) lcolor(white)) bgcolor(white) ///
name(hteHIGH_`var', replace)

graph display hteLOW_`var'
graph save "$paperpath\hteLOW_`var'.png", replace
graph export "$paperpath\hteLOW_`var'.png", replace

graph display hteMID_`var'
graph save "$paperpath\hteMID_`var'.png", replace
graph export "$paperpath\hteMID_`var'.png", replace

graph display hteHIGH_`var'
graph save "$paperpath\hteHIGH_`var'.png", replace
graph export "$paperpath\hteHIGH_`var'.png", replace

}

***************************************************************************************************
***  [Table 4: Heterogeneous Treatment Effect Coaching Conditional on Pre-Coaching Productivity]
***************************************************************************************************

foreach var in $abcd {
	
local verklvar "MonthsToCoachingmin2 MonthsToCoachingplus0 MonthsToCoachingplus1 MonthsToCoachingplus2 MonthsToCoachingplus3"

label variable MonthsToCoachingmin2  "42 -- 22 days before"
label variable MonthsToCoachingplus0 "Day of first coaching"
label variable MonthsToCoachingplus1 " 1 -- 21 days  after" 
label variable MonthsToCoachingplus2 "22 -- 42 days  after" 
label variable MonthsToCoachingplus3 "43 -- 63 days  after" 


esttab hteLOW`var'2 hteMID`var'2 hteHIGH`var'2 using "$paperpath/TABhte`var'.tex", replace f ///
keep(`verklvar') ///
label booktabs b(3) p(3) eqlabels(none) collabels(none) mlabels(none) ///
starlevels($^{*}$ 0.10 $^{**}$ 0.05 $^{***}$ 0.01) ///
cells(b(star  fmt(%9.3f) ) se(par fmt(%9.3f))) ///
stats(N controls driverfe dayfe busdayfe, fmt(%9.0f) labels("Number of trip-level observations" "Controls" "Driver fixed effects" "Day fixed effects" "Bus type $\times$ day fixed effects")) ///
varlabels(_cons Constant) ///
prehead("\tabcolsep=0.25cm" "\begin{tabular}{l*{@M}{C{2.0cm}}}\hline\hline" "Dependent variable: &\multicolumn{1}{c}{\textbf{Lower 25\%}}&\multicolumn{1}{c}{\textbf{Middle 50\%}}&\multicolumn{1}{c}{\textbf{Upper 25\%}} \\" "\cline{2-4}") posthead(\hline) ///
prefoot(\hline) postfoot("\hline" "\end{tabular}")

}

graph close
log close

